function U = singleforce_greensfunc(mdl,obs,mu,nu)
%coded by josh crozier 6/12/2018
%single force greens function
%outputs displacement U
%same input format as disloc3d okada code, except:
    %mdl = [sourcedepth sourceeast sourcenorth eastforce northforce vertforce]
%functions from Segall EQ and Volcano deformation pg 77-78

%forces
f1 = mdl(4);
f2 = mdl(5);
f3 = mdl(6);

%source loc
e1 = mdl(2);
e2 = mdl(3);
e3 = -mdl(1);

%obs locations
x1 = obs(1,:)-e1;
x2 = obs(2,:)-e2;
x3 = obs(3,:);

%source dist from obs
r1 = sqrt(x1.^2 + x2.^2 + (x3-e3).^2);
%imaginary source dist from obs
r2 = sqrt(x1.^2 + x2.^2 + (x3+e3).^2);

%vert greens functions
g13 = x1/(16*pi*mu*(1-nu)).*...
    ((x3-e3)./(r1.^3) +...
    ((3-4*nu)*(x3-e3))./(r2.^3) +...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3)) +...
    (6*x3*e3.*(x3+e3))./(r2.^5));

g23 = x2/(16*pi*mu*(1-nu)).*...
    ((x3-e3)./(r1.^3) +...
    ((3-4*nu)*(x3-e3))./(r2.^3) +...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3)) +...
    (6*x3*e3.*(x3+e3))./(r2.^5));

g33 = 1/(16*pi*mu*(1-nu))*...
    ((3-4*nu)./(r1) +...
    (5-12*nu+8*nu^2)./(r2) +...
    (x3-e3).^2./(r1.^3) +...
    ((3-4*nu)*(x3-e3).^2 - 2*x3*e3)./(r2.^3) +...
    (6*x3*e3.*(x3+e3).^2)./(r2.^5));

%east greens functions
g11 = 1/(16*pi*mu*(1-nu))*...
    ((3-4*nu)./(r1) +...
    1./r2 +...
    (x1.^2)./(r1.^3) + ...
    ((3-4*nu)*x1.^2)./(r2.^3) +...
    (2*(r2.^2 - 3*x1.^2).*x3*e3)./(r2.^5) + ...
    (4*(1-nu)*(1-2*nu)*(r2.^2 - x1.^2 - r2.*(x3+e3)))./(r2.*(r2-x3-e3).^2));

g21 = (x1.*x2)/(16*pi*mu*(1-nu)).*...
    (1./(r1.^3) +...
    (3-4*nu)./(r2.^3) -...
    (6*x3*e3)./(r2.^5) -...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3).^2));

g31 = x1/(16*pi*mu*(1-nu)).*...
    ((x3-e3)./(r1.^3) +...
    ((3-4*nu)*(x3-e3))./(r2.^3) -...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3)) -...
    (6*x3*e3.*(x3 + e3))./(r2.^5));

%north greens functions
g12 = (x2.*x1)/(16*pi*mu*(1-nu)).*...
    (1./(r1.^3) +...
    (3-4*nu)./(r2.^3) -...
    (6*x3*e3)./(r2.^5) -...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3).^2));

g22 = 1/(16*pi*mu*(1-nu))*...
    ((3-4*nu)./(r1) +...
    1./r2 +...
    (x2.^2)./(r1.^3) + ...
    ((3-4*nu)*x2.^2)./(r2.^3) +...
    (2*(r2.^2 - 3*x2.^2).*x3*e3)./(r2.^5) + ...
    (4*(1-nu)*(1-2*nu)*(r2.^2 - x2.^2 - r2.*(x3+e3)))./(r2.*(r2-x3-e3).^2));

g32 = x2/(16*pi*mu*(1-nu)).*...
    ((x3-e3)./(r1.^3) +...
    ((3-4*nu)*(x3-e3))./(r2.^3) -...
    (4*(1-nu)*(1-2*nu))./(r2.*(r2-x3-e3)) -...
    (6*x3*e3.*(x3 + e3))./(r2.^5));

%displacements
U = zeros(size(obs));
U(1,:) = f1*g11 + f2*g12 + f3*g13;
U(2,:) = f1*g21 + f2*g22 + f3*g23;
U(3,:) = f1*g31 + f2*g32 + f3*g33;
end

